Cloud Accretion Simulation: The “Warm Rain” Process#

This notebook simulates the growth of cloud droplets into raindrops through the process of accretion (also known as collision-coalescence). This mechanism is the primary driver of rain formation in “warm clouds” (clouds where temperatures are above freezing), as detailed in the course notes.

The Physics of Accretion#

In a cloud, droplets of different sizes exist simultaneously. Because larger droplets are heavier, they fall faster than smaller droplets. This difference in terminal velocity allows larger drops to overtake and collide with smaller ones in their path.

Before we start to introduce the code, let’s quickly load the modules we will need:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from dataclasses import dataclass
from scipy.spatial.distance import cdist
import math

### from mpl_toolkits.mplot3d import Axes3D

# Enable interactive plots for animation in Jupyter
%matplotlib notebook
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 6
      4 from IPython.display import HTML
      5 from dataclasses import dataclass
----> 6 from scipy.spatial.distance import cdist
      7 import math
      9 ### from mpl_toolkits.mplot3d import Axes3D
     10 
     11 # Enable interactive plots for animation in Jupyter

ModuleNotFoundError: No module named 'scipy'

1. Terminal Velocity#

A droplet’s terminal fall speed (\(w\)) depends on its radius (\(r\)).

We can simplify the relationship by allocating the behaviour to one of three distinct regimes:

  • Stokes’ Law (\(r < 30\mu m\)): \(w \propto r^2\). Viscous forces dominate.

  • Intermediate Regime (\(30\mu m < r < 1mm\)): \(w \propto r\).

  • Newtonian Regime (\(r > 1mm\)): \(w \propto \sqrt{r}\). Aerodynamic drag dominates.

Let’s define these relationships in the following function…

The constants X1, X2, and X3 represent the coefficients for the different fall-speed regimes.

# Physical constants for terminal velocity calculations
X1, X2, X3 = 1.2e8, 8e3, 250. 

def terminal_w(r):
    """Calculation of droplet terminal fallspeed based on radius."""
    # Stokes regime (small drops)
    w = np.where(r < 30.e-6, X1 * r**2, X2 * r)
    # Aerodynamic drag regime (large drops)
    w = np.where(r > 1.e-3, X3 * np.sqrt(r), w)
    return w

Now we need a function for converting radius to drop mass and back again.

This is straightforward, given the density of liquid water \(\rho_l\), the relationship is simply mass=density * volume or

\( M_{drop}=\frac{4}{3} \rho_l \pi r^3\)

rho_l = 1000. # density of liquid water (kg/m^3)

def r2mass(r):
    """Converts radius to mass."""
    return (4.0 * rho_l * np.pi * r**3) / 3.0

def mass2r(m):
    """Converts mass to radius."""
    r = ((3.0 * m) / (4.0 * rho_l * np.pi))**(1./3.)
    return r

2. Collision and Coalescence#

When a large drop overtakes a smaller one, a collision may occur. The probability of these drops sticking together is defined by the Collision Efficiency (\(E\)).

If coalescence occurs:

  1. The masses are combined: \(M_{new} = M_{large} + M_{small}\)

  2. The new radius is calculated from the new mass.

  3. The terminal velocity is updated (the drop accelerates).

  4. The smaller drop is removed from the simulation (“dead”).

The Cloud Class#

This class manages the state of the cloud. It stores the position \((x, y, z)\) and radius \(r\) of every droplet.

The collision method is the computational core:

  1. It calculates the pairwise distance between drops in the X-Y plane.

  2. It checks if drops have “swapped” Z-positions within a timestep (meaning the faster one overtook the slower one).

  3. If they overlap and overtake, a random check against efficiency \(E\) determines if they merge.

Rainfall Calculation#

In the final method of this class, we calculate the rainfall. When a drop passes below the vertical altitude \(z=0\), it is counted as rainfall and removed from the cloud.

class Cloud:
    """Class of a cloud with arrays of x, y, z, r, and w."""
    def __init__(self, xpos, ypos, zpos, radius, collision_E):
        self.x = xpos
        self.y = ypos
        self.z = zpos
        self.r = radius
        self.dead = np.zeros(len(radius), dtype=bool)
        self.rain = 0.0
        self.collision_E = collision_E
        
        # Initialize vertical velocity based on radius
        self.w = terminal_w(radius)   
        self.znew = np.copy(self.z)

    def collision(self):
        ndrop = len(self.r)

        # Calculate pairwise sum of radii (collision threshold)
        # We use a trick to avoid a double loop for speed, though it is memory intensive for large N
        rsq = np.tile(self.r, (ndrop, 1))
        rsum = rsq + np.transpose(rsq) 

        # 1. Check Horizontal Proximity
        coords = np.column_stack((self.x, self.y))    
        xydist = cdist(coords, coords)
        
        # 2. Check Vertical Overtake
        # Did they swap z-locations during this timestep?
        zd1 = np.subtract.outer(self.z, self.z)
        zd2 = np.subtract.outer(self.znew, self.znew)
        
        # Overtake condition: Relative vertical distance changed sign
        # Collision condition: XY dist < radius sum AND overtake happened
        ind1, ind2 = np.where((xydist < rsum) & (zd1 * zd2 < 0.0))
        
        # Filter unique pairs and ignore self-collision
        unique = (ind1 < ind2)
        ind1 = ind1[unique]
        ind2 = ind2[unique]

        # Process collisions
        for i1, i2 in zip(ind1, ind2):
            if self.dead[i1] or self.dead[i2]:
                continue
            
            # Stochastic collision efficiency check
            if np.random.uniform(low=0, high=1) > self.collision_E: 
                continue

            # CONSERVATION OF MASS: Merge 2 into 1
            mass1 = r2mass(self.r[i1])
            mass2 = r2mass(self.r[i2])
            
            # Update Drop 1 (The Collector)
            self.r[i1] = mass2r(mass1 + mass2)
            self.w[i1] = terminal_w(self.r[i1])
            self.znew[i1] = min(self.znew[i1], self.znew[i2]) # Take the lower position

            # Update Drop 2 ( The Collected / Dead)
            self.dead[i2] = True
            self.r[i2] = self.w[i2] = 0.0
            self.znew[i2] = self.z[i2] = -1.e6 # Move well below ground

    def raincalc(self, xysize):
        """Check for drops reaching the surface (z < 0)."""
        rain_idx = np.argwhere(self.z < 0.0)
        
        # Calculate mass hitting the ground
        mass_vec = r2mass(self.r[rain_idx])
        rain_amount = np.sum(mass_vec) / xysize
        self.rain += rain_amount
        
        # Remove rain drops from simulation
        self.dead[rain_idx] = True
        self.r[rain_idx] = 0.0 
        self.w[rain_idx] = 0.0

Part 3: Simulation Parameters & Initialization#

Here we set the size of the domain, the time steps, and the initial population of droplets in a python configuration class, a nice way to handle simulation options as you can also specify the option type.

Experiment Parameters:

  • bimodal: If True, we insert a few large drops manually to kickstart accretion. If False, we use a Gamma distribution (typical for real clouds).

  • ccn: Cloud Condensation Nuclei concentration. Higher values = more, smaller drops (polluted air).

  • collision_E: Efficiency. 1.0 means every collision results in a merger.

# ===========================
# CONFIGURATION CLASS
# ============================
@dataclass
class SimOptions:
    """Holds all parameters for the simulation."""
    # Domain & Time
    xsize: float = 0.0001
    ysize: float = 0.0001
    zsize: float = 2000.0
    dt: float = 60.0         
    lentime: int = 180       
    
    # Physics
    zcloud1: float = 1500.0  
    zcloud2: float = 2000.0  
    collision_E: float = 0.3 
    ccn: float = 200.0       
    rsmall: float = 10.e-6   
    rsigma: float = 5.e-6    
    
    # Distributions
    bimodal: bool = False
    rlarge: float = 25.e-6
    ratio_large: float = 0.05
    
    # Visualization Optimization
    slice_step: int = 2     # Draw only every Nth drop
    steps_per_frame: int = 2 # Run physics N times per video frame

    @property
    def xysize(self):
        return self.xsize * self.ysize
# ==========================================
# INITIALIZATION ROUTINE
# ==========================================

def initialize_simulation(opts: SimOptions):
    """Sets up the drops and creates the Cloud object."""
    
    ccn_m3 = opts.ccn * 1.e6
    ndrop = int(ccn_m3 * opts.xysize * (opts.zcloud2 - opts.zcloud1))
    
    print(f"Initializing Cloud: {ndrop} drops (CCN={opts.ccn})")
    
    if opts.bimodal:
        nlarge = int(np.ceil(opts.ratio_large * ndrop))
        dropr = np.full(ndrop, fill_value=opts.rsmall)
        dropr[-nlarge:] = opts.rlarge
        print(f"- Mode: Bimodal ({nlarge} large drops added)")
    else:
        gshape = (opts.rsmall / opts.rsigma)**2
        gscale = opts.rsmall / gshape
        dropr = np.random.gamma(gshape, scale=gscale, size=ndrop)
        print(f"- Mode: Gamma Distribution")

    # Stats for plotting
    stats = {
        'rmean': np.mean(dropr) * 1e6,
        'rmax': np.max(dropr) * 1e6,
        'rainavail': np.sum(r2mass(dropr)) / opts.xysize,
        'ndrop': ndrop
    }

    # Positions
    dropx = np.random.uniform(0, opts.xsize, ndrop)
    dropy = np.random.uniform(0, opts.ysize, ndrop)
    dropz = np.random.uniform(opts.zcloud1, opts.zcloud2, ndrop)

    cloud = Cloud(dropx, dropy, dropz, dropr, opts.collision_E)
    
    return cloud, stats
def run_visual_simulation(cloud, opts: SimOptions, stats):
    """Runs the loop and generates the animation."""
    
    print("Starting Physics & Animation Loop...")
    %matplotlib inline
    plt.rcParams['animation.embed_limit'] = 50.0

    fig = plt.figure(figsize=(8, 8))
    
    # Setup Axes
    ax0 = fig.add_subplot(211, projection='3d')
    ax0.set_xlim(0, opts.xsize * 1000.)
    ax0.set_ylim(0, opts.ysize * 1000.)
    ax0.set_zlim(0, opts.zsize)
    ax0.set_title(f'Cloud Sim (CCN={opts.ccn})')
    ax0.set_xlabel("x (mm)"); ax0.set_ylabel("y (mm)"); ax0.set_zlabel("z (m)")

    ax1 = fig.add_subplot(212)
    ax1.set_xlim(0, opts.lentime)
    ax1.set_ylim(0, math.ceil(stats['rainavail']))
    ax1.set_xlabel("Time (mins)"); ax1.set_ylabel("Rainfall (mm)")

    # Init Plots
    line, = ax1.plot([], [], 'r.', ms=5)
    ax1.plot([0, opts.lentime], [stats['rainavail'], stats['rainavail']], 'k--', label="Max Rain")
    ax1.legend(loc='upper right')
    
    info = f"Mean R: {stats['rmean']:.1f}um\nMax R: {stats['rmax']:.1f}um\nE: {opts.collision_E}"
    ax1.text(2, stats['rainavail']/2, info, ha='left')

    sc3d = ax0.scatter(cloud.x[::opts.slice_step]*1000., 
                       cloud.y[::opts.slice_step]*1000., 
                       cloud.z[::opts.slice_step], 
                       c='blue', s=cloud.r[::opts.slice_step]*1.e6)

    # State containers
    history = {'times': [0], 'rains': [0.0]}
    sim_state = {'time': 0.0}

    def animate(frame):
        # Physics Sub-stepping
        for _ in range(opts.steps_per_frame):
            cloud.znew = np.copy(cloud.z) - opts.dt * cloud.w
            cloud.collision()
            cloud.z = np.copy(cloud.znew)
            cloud.raincalc(opts.xysize)
            sim_state['time'] += opts.dt

        # Update Graphics
        history['rains'].append(cloud.rain)
        history['times'].append(sim_state['time'] / 60.)
        
        sc3d._offsets3d = (cloud.x[::opts.slice_step]*1000., 
                           cloud.y[::opts.slice_step]*1000., 
                           cloud.z[::opts.slice_step])
        sc3d._sizes = (cloud.r[::opts.slice_step]*1.e6)
        
        live = int(stats['ndrop'] - np.sum(cloud.dead))
        ax0.set_title(f"Time={int(sim_state['time']/60)}m, Live Drops={live}")
        line.set_data(history['times'], history['rains'])
        return sc3d, line

    frames = int(opts.lentime * 60 / (opts.dt * opts.steps_per_frame))
    ani = animation.FuncAnimation(fig, animate, frames=frames, interval=100, blit=False)
    
    plt.close() # Prevent static plot
    return HTML(ani.to_jshtml())
def run_cloud_experiment(**kwargs):
    """
    Main entry point.
    Usage: run_cloud_experiment(ccn=500, bimodal=True, zcloud1=1000)
    """
    # 1. Setup Options (Override defaults with user arguments)
    opts = SimOptions(**kwargs)

    # --- PRINT CONFIGURATION REPORT ---
    print("="*40)
    print(f"{'CLOUD SIMULATION CONFIGURATION':^40}")
    print("="*40)
    for key, value in opts.__dict__.items():
        print(f"{key:<20} : {value}")
    print("-" * 40)

    # 2. Initialize
    cloud, stats = initialize_simulation(opts)
    
    # 3. Run & Animate
    video = run_visual_simulation(cloud, opts, stats)
    
    return video

Part 4: Run Simulation & Animate#

The code below creates a 3D animation.

  • Top Panel: The 3D view of the cloud. The size of the blue dots represents the droplet radius.

  • Bottom Panel: The accumulated rainfall over time.

Watch how the “rain” line stays flat initially (while drops are small and falling slowly), then spikes upwards as larger drops form via accretion and fall out quickly.

run_cloud_experiment()
========================================
     CLOUD SIMULATION CONFIGURATION     
========================================
xsize                : 0.0001
ysize                : 0.0001
zsize                : 2000.0
dt                   : 60.0
lentime              : 180
zcloud1              : 1500.0
zcloud2              : 2000.0
collision_E          : 0.3
ccn                  : 200.0
rsmall               : 1e-05
rsigma               : 5e-06
bimodal              : False
rlarge               : 2.5e-05
ratio_large          : 0.05
slice_step           : 2
steps_per_frame      : 2
----------------------------------------
Initializing Cloud: 1000 drops (CCN=200.0)
- Mode: Gamma Distribution
Starting Physics & Animation Loop...

It is now super straightforward to run sensitivity tests…

# An example sensitivity test....
run_cloud_experiment(ccn=600, bimodal=True, collision_E=0.5)
========================================
     CLOUD SIMULATION CONFIGURATION     
========================================
xsize                : 0.0001
ysize                : 0.0001
zsize                : 2000.0
dt                   : 60.0
lentime              : 180
zcloud1              : 1500.0
zcloud2              : 2000.0
collision_E          : 0.5
ccn                  : 600
rsmall               : 1e-05
rsigma               : 5e-06
bimodal              : True
rlarge               : 2.5e-05
ratio_large          : 0.05
slice_step           : 2
steps_per_frame      : 2
----------------------------------------
Initializing Cloud: 3000 drops (CCN=600)
- Mode: Bimodal (150 large drops added)
Starting Physics & Animation Loop...

Exercises#

To better understand the sensitivity of cloud physics, try the following exercises by modifying the parameters in Part 3.

Exercise 1: The Effect of Cloud Depth#

In the default code, the cloud is relatively shallow (1500m to 2000m).

  1. Hypothesis: What happens if we make the cloud deeper (e.g., extending from 500m to 2000m)? Will rain start sooner or later? Will the raindrops be larger or smaller?

  2. Experiment: Run the code below, which initializes a deeper cloud layer.

Why this matters: Deep clouds allow droplets to fall for a longer duration through the “collection zone,” giving them more time to sweep up smaller droplets before hitting the ground.

Exercise 2: Pollution (CCN Count)#

Change the ccn parameter in Part 3.

  • Low CCN (Clean Air): Set ccn = 50. There are fewer drops, but they share the same available liquid water, so they start larger.

  • High CCN (Polluted Air): Set ccn = 1000. There are many drops, but they are tiny.

  • Question: Does pollution suppress rainfall? Compare the time it takes for the first rain to appear.